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The linewidth of an atom laser is limited by density fluctuations in the Bose-Einstein condensate 
(BEC) from which the atom laser beam is outcoupled. In this paper we show that a stable spa- 
tial mode for an interacting BEC can be generated using a realistic control scheme that includes 
the effects of the measurement backaction. This model extends the feedback theory, based on a 
phase- contrast imaging setup, presented in [T. In particular, it is applicable to a BEC with large 
interatomic interactions and solves the problem of inadequacy of the mean-field (coherent state) 
approximation by utilising a fixed number state approximation. Our numerical analysis shows the 
control to be more effective for a condensate with a large nonlinearity. 
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I. INTRODUCTION 

In recent years there has been interest in utilising 
Bose-Einstein condensates (BEC) and atom lasers for 
precision metrology [2HB]. In particular, the coherence 
properties of such systems make them ideal for perform- 
ing atomic inter ferometry [7 j. However, research has 
demonstrated that the transverse and longitudinal spa- 
tial modes of a BEC exhibit complicated multimode be- 
haviour [HI [9j. This introduces noise into the system, 
which in turn reduces the precision of atom interferomet- 
ric measurements. In a previous paper we theoretically 
demonstrated that a feedback control scheme utilising a 
phase-contrast type measurement could be used to gen- 
erate a stable spatial mode for a BEC possessing negligi- 
ble interatomic interactions pQ. Importantly, the effects 
of the measurement backaction on the system were in- 
cluded in this model. However, many BEC experiments 
work with condensates that have strong interatomic in- 
teractions. While interactions can be removed in some 
systems [10], this adds an additional layer of complex- 
ity to an experiment. Furthermore, atomic interactions 
can be responsible for some useful phenomena, such as 
four- wave mixing [IT] and the generation of nonclassical 
states [I2j [13J. Semiclassical calculations also indicate 
that nonlinear interactions are necessary for the stabil- 
ity of continuously pumped atom lasers [14] [15] . In this 
paper we further develop the theory presented in [1] to 
show that feedback control can be used to generate a 
stable spatial mode for an interacting BEC. 

It is experimentally possible to create condensates with 
negligible atomic interactions. This can be done by us- 
ing a dilute atomic sample or via a Feshbach resonance 
[16]. The dynamics of a noninteracting trapped BEC is 
very similar to a trapped single atom. Thus work done 
in controlling a single atom is applicable. Doherty and 
Jacobs [17] showed that feedback control could be used 
to stabilise an atom that had its position continuously 
monitored. This was done by solving the optimal control 
problem for an initial Gaussian state. Such a position 



measurement could be engineered by magnetically trap- 
ping the atom in an optical cavity. Wilson et al ex- 
panded upon this work by showing that the stochastic 
master equation (SME) for the model could be solved, 
and thus the atom could be cooled from any arbitrary 
state [18] . However, continuously measuring the atom's 
position required that the atom be trapped in a region 
small compared to the wavelength of light within the cav- 
ity. This condition is not met by a modestly sized BEC 
trapped in an optical cavity. We addressed this issue in 
a previous paper [1 by deriving a control scheme for a 
noninteracting BEC based upon phase-contrast imaging, 
a nondestructive density measurement that has already 
been utilised in experiments [I9j [20] . It was shown that 
in the single atom limit a robust feedback control, based 
upon semiclassical work performed by Haine et al. [21] 
(see also [22 ), would drive the atom towards a stable 
spatial mode close to the ground state energy. 

However, large nonlinearities associated with atomic 
interactions arise naturally in typical BEC experiments 
(for instance [23-25 ). If one wants to design and build 
an atom laser for use in precision metrology, such non- 
linearities cause a number of theoretical and practical 
challenges. It has been demonstrated that in an atom 
laser, interatomic interactions cause number fluctuations 
to couple to energy fluctuations, which broadens the out- 
put beam linewidth [26j[27]. Furthermore, pumping an 
atom laser excites the spatial modes of the lasing mode 
[T4| IT5| [28] . Naively, it may seem that removing inter- 
atomic interactions from the system could be advanta- 
geous. There are, however, a number of reasons why 
it would be preferable to be able to control a conden- 
sate with high interatomic interactions. From a practical 
standpoint, the creation of a noninteracting BEC creates 
an additional layer of experimental complexity. Fesh- 
bach resonances require precise control of the absolute 
magnetic field, and therefore preclude magnetic trapping 
of the condensate. Outcoupling atom lasers is harder 
in optical traps, as optical traps are typically far less 
state- selective. More importantly, theoretical modelling 
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predicts that a continuously pumped atom laser is only 
stable in the regime of high atomic interactions [T4} [T5]. 
More generally, there are situations where the presence 
of interactions results in interesting phenomena worth 
studying for their own sake. Proposals to generate non- 
classical states in atom laser beams [I2j [13], four- wave 
mixing experiments [IT] and the Bosenova experiment of 
Donley et. al [29] are a few examples. 

Some work by Wiseman and Thomsen has shown that 
for a single mode atom laser, feedback control can be used 
to reduce the phase diffusion caused by large atomic in- 
teractions [26l[27]. More recently, Yanagisawa and James 
have proposed using coherent control to directly cancel 
the effects of the phase diffusion [30;. However, such 
schemes do not address the noise associated with insta- 
bility in the BEC spatial mode, which is often the dom- 
inant effect. In this paper, we show theoretically that 
the control setup considered in [1 can be used to drive 
an interacting BEC to a steady state close to the ground 
state energy. 

The structure of the paper is as follows. In Sec. [II] 
we present our full-field model of the system, measure- 
ment apparatus and feedback, and the associated SME 
for the quantum filter. A derivation of this SME can be 
found in the appendix of [1 . In Sec. Ill we simplify the 



quantum filter by making a semiclassical approximation. 
More precisely, we assume that the state vector is always 
in a Fock state of fixed total number (the Hartree ap- 
proximation). Note that the mean- field approximation, 
that has been so successful in BEC theory, is unsuit- 
able in this system since the measurement projects the 
BEC towards a number state. Using numerical simula- 
tions of the order parameter under this approximation, 
we demonstrate that our control scheme does give cooling 



to a steady state. In Sec. IV we argue that on a timescale 
short compared with the time required to reach steady 
state, the semiclassical wavefunction is projected to a 
Gaussian function. Approximating the state thus, we 
perform a numerical analysis on this model. These sim- 
ulations demonstrate that (a) there is an optimal choice 
for feedback, (b) the average steady-state energy scales 
with the measurement strength, and (c) increasing the 
interatomic interaction strength cools the BEC closer to 
the ground state energy. 



II. FULL-FIELD QUANTUM MODEL 

The system under analysis is a BEC magnetically con- 
fined in an harmonic trap (of frequency uj x in the x direc- 
tion and uj±_ in the y and z directions, with uj x <^ujj_) and 
illuminated with off-resonant laser light directed along 
the z-axis (see Fig. [I]). The condensate's density profile 
along the x-axis is obtained from homodyne detection of 
the light after it has interacted with the atoms. The con- 
trol is performed by enacting feedback via adjustments of 
the trapping potential. In [1] we presented a mathemat- 
ical model for this control setup, based on a system-bath 
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FIG. 1. (Colour online) Control setup modelled in this pa- 
per. A BEC confined in an harmonic trap is illuminated with 
an off-resonant coherent light field. A measurement of the 
condensate's density profile along the x-axis is performed us- 
ing homodyne detection in the phase quadrature. Informa- 
tion obtained from this measurement is used to construct an 
estimate of the quantum state of the BEC (eq. |l])). This es- 
timate is used to perform real-time feedback on the BEC via 
the adjustment of the magnetic trapping potential (eq. Q). 



coupling between the trapped BEC and the electromag- 
netic field. We then derived from this model the following 
conditional master equation: 

dp c = — i[H, p c ]dt + a J dxV[M(x)]p c dt 

+ yfe j dxH[M(x)]p c dW(x,t). (1) 

For convenience we have expressed position and time in 
dimensionless harmonic oscillator units, where distance 
is in units of xo = ^h/muj x and time is in units of cj" 1 . 
Here m is the mass of the atomic species. The conditional 
density operator p c is the best estimate (in the least- 
squares sense) of the quantum state of the BEC [31]. 
Now let us consider each individual term of eq. 0. The 
first term describes the unitary dynamics of the system, 
which has Hamiltonian 



H 



U d 



dxip^ (x)H a (x)ip(x)-\-^- J dx ffi (x)ffi (x)$(x)i[)(x) , 

(2) 

where x/j(x) is the field operator that annihilates an atom 
in the ground state at position x (the detuning A is suf- 
ficiently large that we have adiabatically eliminated the 
excited state). The field operators obey the commuta- 
tion relation [^(x), ^ (x 1 )} = 8{x — x'). The first term 
in eq. ([2| is the single atom Hamiltonian, containing the 
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kinetic energy and the trapping potential: 



H a (x) = 



1 d 2 



H, 



control 



(x). 



(3) 



^control (x) is the single particle control Hamiltonian, 
which is caused by adjustments in the x-direction mag- 
netic trapping potential. We consider feedback of the 
form 



control 



(x) = Ci (ft) X + C2 (xp + px) X 2 



(4) 



where c\ and C2 are constants that determine the strength 
of the feedback. The first term in eq. Q models adjust- 
ments to the magnetic trap minimum, while the second 
term represents adjustments made to the 'tightness' (that 
is, the gradient) of the harmonic potential. Thus feed- 
back proportional to x and x 2 will control the 'sloshing' 
and 'breathing' modes of the condensate, respectively. 
It should also be noted that since this is a many body 
system, the operators x and p take the form 



dx ffl (x)x ip(x) 
dxip\x) ) 1 P( X )- 



(5) 
(6) 



The second term in eq. Q models the energy due to 
collisions between the atoms. In gases of ultracold al- 
kali atoms the range of the scattering potential is much 
less than the average spacing between atoms. Thus 
the scattering potential can be adequately modelled by 
a hard-sphere contact potential, as has been done in 
eq. ([2]). Furthermore, at low energies s-wave scatter- 
ing dominates. Hence the interatomic interactions are 
determined by an effective 3D interaction parameter 
U3 = 4:7rh 2 a s /m, where a s is the s-wave scattering length 
[16]. In eq. pi) the dimensionless ID interaction strength 
Ud — (Ui/xo)/huj x is the effective interaction parameter, 
where U\ is roughly Us divided by the transverse area 
of the condensate [32 . The atomic interactions were ne- 
glected in our previous paper pQ. 

The second and third terms in eq. ([!]) are due to the 
interaction between the BEC and the light field. The 
strength of this interaction is given by the measurement 
strength parameter 



a 



3 T S p 
4^7 A2' 



(7) 



where r sp is the rate at which a single atom sponta- 
neously emits into the environment, Q is the Rabi fre- 
quency and A is the detuning of the laser. Notice that 
increasing the intensity of the laser (which increases Q) 
or moving the frequency of the laser closer to the atomic 
transition (decreasing A) results in a larger a. Although 
a larger a gathers more information per measurement 
(the third term of eq. Q), it also increases the rate 
of heating of the atomic ensemble (the second term of 
eq. 0). 



The second term in eq. ([!]) features the decoherence 
superoperator 



V[c]p c = cp c c ] - ^{c^pc}, 



(8) 



where c is any arbitrary operator. This term is the de- 
coherence experienced by the condensate at each point x 
due to the measurement 



M(x) 



dx' T(x - x')^ (x')^(x'), 



(9) 



whence 



T(x) 



2iri]_ 



J dk VT^e 



j(k) = exp 



1 ( V± 

T] 2 



(10) 

(ii) 



Here we have defined the Lamb-Dicke parameter 77 = 
koxo and r]± = koR±, where R± is the length of the 
condensate in the y and z-directions and ko = 2tt/X is 
the wave number of the incoming laser of wavelength A. 
The expression for j(k) (eq. (11)) is only applicable in 
the limit R z ^> A (i.e. r]± ^> 27r), which is the limit 
where photons interacting with the BEC are predomi- 
nantly scattered in the forward z-direction. One can see 
from the measurement operator M(x) that the interac- 
tion between the condensate and light field results in a 
measurement of the number at position x (ffl (x)i/j(x)), 
blurred by the function T(x). Indeed, the width of T(x) 
gives the resolution length scale of the measurement. 
Thus the second term in eq. represents the decoher- 
ence due to the measurement backaction. 

The third term in eq. 0, called the 'innovations', rep- 
resents the new information gathered via the measure- 
ment process. From another perspective, one can think 
of the innovations term as the measurement signal ob- 
tained after homodyne detection. The new information 
obtained about the condensate at each point x from the 
measurement M(x) is encoded in the superoperator 

H[M(x)]p = M(x)p + pM\x) - Tr[(M(x) + M\x))p]p. 

(12) 

The homodyne measurement signal is corrupted by 
quantum noise due to random wave function collapse. 
This noise, within a limited bandwidth, is Gaussian 
white noise and is modelled using the Wiener increment 
dW(x,t). It satisfies dW(x, t)dW(x', t) = 5(x - x')db. 
The Wiener increment vanishes when we take the en- 
semble average - i.e. E[dW(x, t)] = 0. Furthermore, for 
any physical operator c, E[cdW(x,t)} = E[c]E[dW(x,i)]. 
Thus if we have no control (set c\ = C2 = 0) and take the 
ensemble average of SME (HI) , then we recover the master 
equation for p = E[p c ]. 
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III. SEMICLASSICAL MODEL: THE HARTREE 
APPROXIMATION 

Eq. ([!]) contains the full quantum dynamics of the 
BEC. Unfortunately, it is impossible to obtain an ana- 
lytic solution to this equation. Moreover, it is also unfea- 
sible to obtain a solution via numeric integration due to 
the extremely high dimensionality of the quantum field. 
We must therefore make an approximation. In typical 
BEC experiments many of the quantum correlations are 
unimportant, and can be neglected [33 . Indeed, in many 
experiments only lower order moments, such as average 
density, are measured. In such cases, the relevant dynam- 
ics of the condensate can be adequately modelled with 
a semiclassical 'mean-field' wavefunction. In the BEC 
literature, a particularly successful semiclassical approx- 
imation is to assume the BEC is always in a specific co- 
herent state. However, as briefly mentioned in pQ, such 
an approximation is inappropriate for the above control 
apparatus. Primarily, this is due to the type of measure- 
ment we are performing on the BEC. A continuous weak 
measurement of the form ([9| will, over time, project the 
BEC to a fixed global number state. This is not in good 
agreement with a coherent state, which is a Poissonian 
distribution of number states. 

Given these considerations, it is more reasonable to use 
a semiclassical approximation that assumes (a) the BEC 
is always in a number state of fixed total number, TV, 
and (b) all TV atoms occupy the same mode. That is, we 
make the Hartree approximation, where we assume that 



the state can always be written in the form 

|¥) = |JV,0,0,...>, (13) 

in some (possibly time-dependent) single particle basis. 
By writing i/j(x) = ^ n Xn(x)K, where b n are the creation 
and annihilation operators for the above number basis, 
we can see that 

(ft{x)$(x')) = <*|^(aO^(aO|*) =N X *(x)x(x') 

(14) 

and 

(y\x)^(x')^{x)^(x')) = N(N - l)| X (z)| 2 |xOr')| 2 - 

(15) 

x(x) = Xo( x ) is the order parameter for the mode con- 
taining N particles. 

In order to derive an equation of motion for x( x )i we 
cannot simply compute 

d(y(x)) = Tr{$(x)dp c } (16) 

since (i/j(x)) = 0. Instead we must consider the one- 
body density matrix n^ 1 \x,x / ) = (ft ' (x)ip(x')) , which 
has non-trivial evolution: 

d ^ f (x)VV)) = Tr \j)\x)^(x')dp c } . (17) 

Substituting SME ([!]) into eq. ( [l7| ) and performing some 
straightforward operator algebra yields [34] 



dnW(x,x') = -iN {x*(x)H a (x') X (x') - X (x')H a (x) X * (x) + U d (N - 1) (\ X (x')\ 2 - \ X (x)\ 2 ) X *(x) X (x')} dt 
+ N { a f dkj(k) (<><*(*-*')_ dt + ^J ^VtW [(e <fcx - (e ik ^dW(k,t) 

+ ( e -**' - (e- ifc <->)) dW*(k,t)] } X *(x) X (x'), (18) 



,»fe(.) 



dx X *(x) e ikx X (x) 



(19) 



and 



where lations: 

dW\k,t)dW{k',t) =5(k-k')dt (21) 
dW(k, t)dW{k\ t) = 5{k + k')dt. (22) 

We can express dW(k,t) in terms of another Wiener in- 
crement dY(k, t). Specifically, 

dW(k,t) = -(i-l)[dY(k,t) + idY{-k,t)]. (23) 
By applying the Ito product rule to eq. ( [14] ) we obtain 

dnW{x, x f ) = N [x(x f )d X *(x) + x*i?)dx(x') 

__ +dx*(x)d X (x')}. (24) 

dW(k, t) is the Fourier transform of the Wiener incre- 
ment. It is complex- valued, and has the following corre- An equation for dx(x) that satisfies eqs (18) and (24) is 



dW(k,t) 



dke~ lkx dW(x,t). 



(20) 



5 



ik(-) 



d X (x) = !^-iH(x)dt - | J dkj(k) (l - 2e~ ikx 

J ^v^(e- lb -(e-^))dr(M)jxW, 



dt 



(25) 



where 



H(x) 



i d 2 i 

~2~dx^ + 2"' 
+ ci (p) a; - 



2 + (N 
c 2 (xp - 



-i)^lxWI 2 

- px) x 2 



(26) 



is the system Hamiltonian for the BEC. 



Note that the decoherence and innovations in eq. ( 25 ) 



are independent of the total number of atoms TV. In fact, 



the nonlinear term in the system Hamiltonian (26) is the 
only term that depends on N. Furthermore, if we set 
N = 1 then we recover the single atom limit of SME 
(cf. eq. (29) in [1 ), as we would expect. 



A. Simulation of eq. (25) 



The primary aim of our control scheme is to drive the 
BEC towards a stable spatial mode - that is, a steady 
state. We would also like this mode to be (a) close to 
the ground state energy, and (b) obtainable in an experi- 
mentally reasonable period of time. How well the control 
scheme satisfies these three criteria is best judged by ex- 
amining the average energy of the BEC: 

E[E} = ±E[(p 2 ) + (x 2 )+u(\ X \ 2 )], (27) 

where u = (N — l)Ud is the effective interaction strength. 
The energy ([27]) was calculated by finding a numerical so- 
lution to the stochastic Schrodinger equation (SSE) for 
x(x), namely eq. (25). The numerical integration was 
performed with the open source software package xmds2, 
which is a new version of the xmds package [35]. An ex- 



ample of the typical dynamics revealed by solving eq. ( 25 ) 
is shown in Fig. |2j As one can see, a higher energy ini- 
tial state can be cooled to a steady-state of lower average 
energy. As we expect, in the limit of small u, the cool- 
ing rate and the final steady-state energy follow similar 
trends to those in the single atom limit. Our previous 
work in [Tj details these results, including the approxi- 
mate scaling of the final energy as arf . 



IV. THE GAUSSIAN ASSUMPTION 

As stated in the previous section, numerical solutions 
to the SSE (25) indicate that in the limit of small u and 



V± ^ V the control proportional to x cools a highly ex- 
cited state to a steady state. However, we are interested 
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FIG. 2. (Colour online) Plot of the energy per particle of the 
BEC as a function of time for 100 trajectories. The solid line 
denotes the average energy, whilst the dashed lines around the 
solid line indicate the standard error. The energy is in units 
of Eg w 1. 1735 fouj x , which is the ground state energy of the 
condensate. This energy, which was calculated numerically, is 
marked on the plot by the flat dashed line. The parameters 
chosen for this simulation are a = 1.0, r\ — 4.0, rj± — 20.0, 
c\ — 2.0, C2 = 0.0 and u = 4.0. The initial condition was 
a normalised Gaussian wavefunction offset at x = 5 with a 
width a = Ea- 



rn controlling a strongly interacting condensate where 
u is large. Furthermore, in a realistic BEC experiment 
r]± ~ f]. Unfortunately this regime can only be simulated 
on short timescales (compared with the time required 
to reach steady state) due to the current limitations of 
numerical algorithms for stochastic differential equations 
(SDEs) and computational power. Fortunately, there is 
good reason to believe that we can make a further ap- 
proximation to this system and still obtain insightful nu- 
merical results. Fig. [3] shows that over a short period of 
time an initial off-centred Thomas Fermi wavefunction, 
in the large r] ~ r]± limit where both the x and x 2 con- 
trols are on, is driven towards a state that is Gaussian. 

This observation motivates us to assume that the order 
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FIG. 3. (Colour online) Simulation of the ensemble average 
of \x(x,t)\ 2 at times t — (red dashed line) and t — IOcj" 1 
(blue solid line). The state at t = is the ground state 
Thomas-Fermi wavefunction offset to position x = 10 xo- The 
density profile at t — 10 u^ 1 is certainly Gaussian, since a 
least-squares fit of a Gaussian function gives R 2 = 0.995. 
Parameters for the simulation are a = 1.0, r\ — 6.0, rj± = 
10.0 and u = 60.0. The average is over 800 paths. Curves 
indicating the standard error are sufficiently small that they 
have been omitted for the sake of clarity. 



parameter x( x ) is always of the form 



x(M) 



exp 



(iV xp - l)(x - (x)) 2 



■ i (p) X 



(28) 

where V xx is twice the variance in x and V xp is the sym- 
metrised covariance. That is 



Va, 



(<* 2 > - (x) 2 ) 



V xp = (xp + px) -2 (x) (p) . 



(29) 
(30) 



By making this approximation, we have assumed that 
the important dynamics of the system only depend upon 
four variables. The four coupled Ito stochastic equations 
of motion for these variables are (see Appendix [A] for 
details): 



d (x) = (p) dt + y/pV^dW^t) 



(31) 



d(p) = - ((x) + ci (p) + 2c 2 (2 (x) (p) + V xp ) (x)) dt 



(32) 



dV xx = 2 (y xp - PV^I 2 ) dt + ^/3pV^ 4 dW 2 (t) (33) 



2c 2 (2 (x) (p) + V xp ) V x: 




dt 



(34) 



where dW\ (t) and dW^ (t) are independent Wiener incre- 
ments (such that dW\(t)dW2(t) = 0) and j3 = a/rjrj± is 
the effective measurement strength. 



Numerical Results 



As outlined previously in Sec. Ill A[ we determine that 
our control has driven the BEC to a steady state when the 
average energy (eq. (27)) reaches a steady state. Under 
the Gaussian approximation the average energy takes the 
form 



E[E] 



1 



E 



(pY 



xp 



2V r , 



{xY 



i 



^kV~ x , 



(35) 

We numerically solved eqs (31 )-(34) and output the aver- 
age energy (35) using the xmds2 package. In particular, 



we studied the effects of control on the BEC due to the 
four free parameters: the feedback strengths c\ and C2, 
the effective measurement strength fj and the nonlinear 
interaction strength u. The results of this analysis are 
outlined below. 

The first result is that for each /3 and u there exist 
optimal values for the feedback strengths c\ and C2 that 
minimise the average steady-state energy, and the time 
taken to attain the steady state. Figs[4^ and[4]3 show that 
as one varies each individual feedback strength, the aver- 
age steady-state energy goes through a minimum. This is 
most easily understood by noting that the 'sloshing' (ci) 
and 'breathing' (02) controls dampen the (x) and (x 2 ) 
modes of the condensate, respectively. As illustrated in 
Fig.|4j3, the choice of feedback strength gives three differ- 
ent regimes of feedback control: underdamped, critically 
damped and overdamped. Optimal control occurs when 
the feedback strengths are chosen to give critical damping 
of the modes (x) and (# 2 ), as this cools the condensate 
to a minimum energy steady state (for a given f3 and u) 
in the minimum amount of time. For the remainder of 
this analysis, we have chosen c\ and C2 close to optimal. 

Fig. [5] shows the effect on the average steady-state en- 
ergy as the effective measurement strength j3 is varied. 
The trend indicates that a larger j3 results in a higher 
energy steady state. This makes sense, as a larger mea- 
surement strength means that there is a greater measure- 
ment backaction on the condensate. An increased back- 
action increases that rate at which energy is transferred 
from the light field to the atoms in the condensate. Thus, 
while a steady state is still attained, it is of a higher en- 
ergy. The general scaling of the relationship between f3 
and the average steady-state energy is difficult to ascer- 
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FIG. 4. (Colour online) Plots illustrating the existence of optimal feedback parameters for fixed f3 = 0.04 and u = 8.0. (a) 
Numerical solutions to the average steady-state energy as a function of c\ for fixed c 2 = 0.32. The minimum of E[E] ^ 1.817^ 
occurs at c\ — 1.3, where is the ground state energy of the BEC. (b) Numerical solutions to the average steady-state energy 
as a function of c 2 for fixed c\ — 1.30. The minimum of K[E] ~ 1.817^ occurs at ci = 0.39. (c) Numerical simulations of the 
(x 2 ) moment of x( x ) over time £, for fixed c\ — 1.30. Depending on the choice of c 2 this mode is either i. underdamped, ii. 
critically damped or iii. overdamped. The optimal choice is a feedback strength that gives critical damping, as this cools the 
BEC to the lowest steady-state energy in the shortest interval of time. Similar regimes can be shown to exist for (x) by varying 
c\. In (a) and (b) the standard error in the average steady-state energy is smaller than the width of the points that indicate 
the mean, whilst in (c) the mean and standard error in (x 2 ) are given by the solid and dashed lines, respectively. 



tain numerically, since (as is shown below) it is dependent 
upon the interaction strength u. 

The final result concerns how the size of the inter- 
atomic interaction strength affects the control. Fig. [6] 
shows how the percentage difference between the average 
steady-state energy and the ground state energy changes 
as u is varied. The general trend is that as u increases the 
energy difference decreases, and thus the steady state is 
closer to the ground state energy. Specifically, the steep- 
est decrease in the energy occurs when the interaction 
strength is increased from zero to a medium strength 
(u = 16-32). For instance, for f3 = 0.04 the energy dif- 
ference decreases by about 6% when u changes from 
to 16. In contrast, increasing u from 128 to 512 only 
decreases the energy difference by 0.55%. The increased 
effectiveness of the control at removing energy from in- 
teracting condensates, compared to noninteracting con- 
densates, can be explained by studying which modes are 



affected by the feedback. Recall that the 'sloshing' and 
'breathing' controls remove energy from the (x) and (x 2 ) 
modes of the BEC, respectively. However, these controls 
do not directly remove energy from higher order modes 
((x 3 ), (^ 4 ), etc.). Hence, in the limit of no atomic inter- 
actions, these higher order modes are unaffected by the 
control and remain excited. However, for a condensate 
with atomic interactions the nonlinearity couples (x) to 
higher order odd modes, and (x 2 ) to higher order even 
modes. Thus, through this coupling, the feedback con- 
trols can dampen these higher order modes. Note that 
the effect of this coupling is largest in the low u regime, 
as is indicated by the steep decline in energy for small 
changes in u (see Fig. [6|. In contrast, while increases 
in the interaction strength for larger values of u still de- 
crease the energy difference, this effectiveness is reduced 
simply because there is a finite amount of energy stored 
in the higher order modes. Stronger nonlinearities may 
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FIG. 5. (Colour online) Plot showing how the average steady- 
state energy compares to the ground state energy as a function 
of the effective measurement strength /3. This demonstrates 
that a larger choice of f3 results in a higher energy steady 
state. Each data point was generated from a numerical simu- 
lation where optimal feedback parameters were chosen. The 
standard error is smaller than the width of each point. 



'squeeze' out additional bits of energy from these modes, 
however the vast majority is removed in the low u cou- 
pling regime. Note also, from Fig.[6j that the sharp drop 
off in the average steady-state energy difference is more 
pronounced for f3 = 0.08 than j3 = 0.04. Indeed, for 
/3 = 0.08 the energy difference decreases by 17.6% when u 
changes from to 16, compared with 6% for that same u 
interval. One could conclude, therefore, that the presence 
of nonlinearities in the condensate are more important 
for effective control in the regimes of larger measurement 
strength. 



B. Remarks on validity of Gaussian assumption 

Many of the calculations presented in this paper are 
based upon the assumption that the state x( x ) evolves 
towards a Gaussian function. However, this can only be 
proved for a Linear-Quadratic-Gaussian (LQG) system 
(see, for example, [U]). The Gaussian assumption re- 
duces the total number of modes that are affected by 
vacuum noise and heating. By making this assumption, 
we have implicitly assumed that most of these channels 
of noise are small compared with the ones that strongly 
affect a Gaussian state, and can thus be neglected. We 
checked that the state approaches a Gaussian by making 
numerical simulations of the complete semiclassical equa- 
tion (cf. Fig. [3|. These numerical checks were performed 
for appropriate physical parameters, and showed that the 
state reached a Gaussian on average after several trap cy- 
cles and stayed Gaussian thereafter. The approximation 
then allowed the simulations to be run over much longer 
timescales, so that the final steady-state energies could 
be determined. 



FIG. 6. (Colour online) Plot showing how the average steady- 
state energy compares to the ground state energy as a function 
of the interaction strength u for (blue dot) f3 = 0.04 and 
(maroon square) /3 = 0.08. One can see that the control 
cools the BEC closer to the ground state for larger interaction 
strengths. Each data point was generated from a numerical 
simulation where optimal feedback parameters were chosen. 
The standard error is smaller than the width of each point. 



V. CONCLUSIONS 

We have presented a model of an interacting BEC 
undergoing feedback control via a continuous dispersive 
imaging measurement. More precisely, we considered 
the filtering equation derived in pQ under the semiclas- 
sical Hartree approximation. This approximation, where 
the BEC is assumed to be in a number state, was used 
in preference to the mean-field approximation since the 
measurement projects the BEC into a single number 
state, rather than a coherent state. Numerical simu- 
lations showed that the mean-field approximation gave 
unphysical results, whereas the Hartree approximation 
gave physically appropriate dynamical behaviour. We 
further refined our semiclassical model by assuming that 
the semiclassical wavefunction x( x ) was always a Gaus- 
sian function. We then provided a numerical analysis 
of this model, where the affects of the free parameters 
on the effectiveness of the control scheme were studied. 
There were three key results from this analysis: 

1. For each u and /?, there exist optimal values for 
the feedback strengths that minimise the average 
steady-state energy and the time taken to attain 
the steady state (cf. Fig. |4|; 

2. The average steady-state energy increases as the 
effective measurement strength j3 increases (cf. 
Fig.§; 

3. The average steady-state energy (relative to the 
ground state energy) decreases with increasing 
atomic interaction strength u (cf. Fig. [6]). 
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The final point is of particular interest, as it indicates 
that the control scheme is more effective for a strongly 
interacting BEC, which is the situation in many BEC 
laboratories. Thus, if one wanted to implement this con- 
trol scheme, there is no need to expend effort and re- 
sources removing the nonlinearities of the BEC. Most 
importantly, our work has shown that this is a viable 
control scheme for reducing the multimode fluctuations 
of a trapped BEC. 
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Appendix A: Derivation of eqs (31) - (34) 



Before we apply the Gaussian ansatz (28) it is neces- 
sary to calculate the equations of motion for conditional 
expectation values. From the Ito product rule we know 
that the dynamical equation for the expectation value of 
x, for example, under the Hartree approximation is 



d (x) — d ( / dx x*(x)%x( x ) 



I 



dxx[dx*(x)x(x)+x*(x)dx(x) + dx*(x)dx(x)} . 

(Al) 



Substituting eq. ( 25 ) into equations of m otion for expec- 
tation values (similarly expressed as eq. (Al)) yields 

d (x) = (p)dt+y/a / dk VrtX) [A(k)dW k + c.c] 



(A2) 



d(p) 

d(x 2 ) 
d( P 2 ) 



d (xp + px) = 2 ((p 2 ) ~(x 2 )- ci (a:) (p) 

-2c 2 (xp + px) (x 2 )) dt + iu ((\x\ 2 px) 

- (px\x\ 2 )) dt + 2v^y dk yffik) (D(k)dW k 
+ (D(-k)-kC x {k))dW* k ) (A6) 



:((x) 2 ) = 2 (x) (p) dt + 4a / dk 1 (k)A(k)A(-k)dt 
+ 2v^ (x) J dk ^/^(k~j [A(k)dW k 



C.C.I 



(A7) 

d ((p) 2 ) = -2 (p) ((x) + ci <p> + 2c 2 (xp + px) (x)) dt 
+ 2^(p) J dk^(k)[B(k)dW k + (B(-k) 

-kC(k))dW* k +aj dkj(k)[4B(k)B(-k) 
-AkB(k)C(k) - k 2 C(k)C(-k)] dt (A8) 
d((x)(p))=((p) 2 -(x) 2 -c 1 (x) (p) 

— 2c2 (xp + px) (x) 2 ^j dt 

+ 2c*y dk-f(k) [A(k)B(-k) + A(-k)B(k) 

-kA(k)C(k)] dt 
+ y^Jdk y/rfe) [<P> {A(k)dW k + c.c} + 
(x) {B(k)dW k + (B(-k) - kC(k))dW* k ^ , 



(A9) 



where 



: - ((x) + ci (p) + 2c 2 (xp + px) (x)) dt 

J dk B(k)dW k + (B(-k) - kC(k))dW] 



A(k) = 


(xe ikx ) - 


(x) (e ikx ) 


A 2 (k) = 


(x 2 e zkx ) - 


- (x 2 ) (e* kx ) 


B(k) = 


(e ikx p) - 


(P) (e ikx ) 


B 2 (k) = 


(e* kx p 2 ) - 


-(P 2 )(e^) 


C(k) = 


(e~ ikx ) 




C x (k) = 


(xe~ ikx ) 




D(k) = 


(e ikx xp) 


- (xp) (e ikx ) 



(A3) 

(xp + px) dt + \fa J dk V7(fe) [A 2 (k)dW k + c.c] 

(A4) 
dt 



After assuming the ansatz ( |28| ) we obtain after some sim- 
plification: 



d (x) = (p) dt 



arj 



(^(xp + px) + 2ci (p) 2 + 2c2 (xp + px) 
dk k 2 'y(k)dt + ^ 



2(2tt 



,1/4 



Vt, 



dx \F\{x) — c.c] dW(x,t) 
(A10) 



dk^Mk)[B 2 (k)dW k d(p) 



+(B 2 (-k) + k 2 C{k) - 2kC x (k))dW k 



((x) + a (p) + 2c 2 (2 (x) (p) + V xp ) (x)) dt 
1 / a 



(A5) 



2(2tt) 1 /4 y V j_ 



dx [(1 + iV xp )F 1 (x) + c.c] dW(x, t) 

(All) 



10 



dV xx = 2V xp dt 



2 a 



2tt V± 

(27T) 1 /4 f^ V 2 



V 2 x g(t)dt 



= / dx \F 2 (x) + c.c.1 dWYx, t) 

2 V T)± T) J 

(A12) 



dV xp = 



VI + 1 



xp 



V xx - 2c 2 (2 (x) (p> + V xp ) V xx dt 



XX 

U 1 



-^-2V2tt— V^K^)^ 

V27T V Via; V± 



v XX 

1 ra~v T . 



2(2tt) 1 /4 y ^ ^ 

where we have defined 
1 



*i(aO 
F 2 {x) 



2tt 
1 



d/c /c exp 
dkk 2 exp 



xdW(x,i), (A13) 



^l^ 4 _ ^KcaA ik(x/rj-(x)) 

1677 4 4 y 

(A14) 

^l^ 4 _ ^^aA ik(x/rj-(x)) 

1677 4 4 y 

(A15) 
(A16) 



ee / dfc A: 2 exp ( - T ^~ ~ 



Note that we have expressed the equations of motion in 
terms of the real, x-space noise dW(x,t). 

In order to proceed further it is necessary to make an 
approximation. We are going to choose para mete rs such 
that the quartic term in the exponent of eqs (A14)-(A16) 
can be neglected. That is 



F 2 (x) 



2tt 
2x/2i 

T7 4 V 3/2 
'/ * XX 



1 r dkkexp(-^^)e^^-^ 



(x-r]{x))e-^-^ x))2/ri2Vx - 
dkk 2 exp 



2tt 
2y/2 



5/2 

XX 



(A17) 



( V 2 V XX -2(x-r) (x)) 2 ) e -(*-v(x)fh*v xx 



1 f " -2 ( ] e i k (,:/,,-(x)) 



I 



g(t) & dkk exp 



2 / ^ v xx 



2tt 



t/3/2 * 



(A18) 



(A19) 



This approximation is valid when the fc-space variable in 



the integrands of eqs (A14)-(A16) satisfy 



fc 2 « 



\/k G [ ^neg? ^neg] 5 



(A20) 



where the integrand in the above integrals is negligible 
(i.e. roughly zero) at k = ±& neg . 



Under this approximation, eqs (A10)-(A13) can be 
written as 

d (x) = (p) dt 

2 5/4 V^ f , ,(x'-rj(x)) -i^M) 



/ 



dx' v ~ ^T —e ^v xx dW(x',t) 

(A21) 



d(p} = - ((x) + a (p) + 2c 2 (2 (x) <p) + K P ) <ar)) dt 

->5/4./w, r , V xp (x' - T] (x)) -li^iMl 
~i — o — ; — ; "«£ " ' 

7T 



= 2V xp dt 



2 b ^^ r 

2a vH 2 dt 

/• 



yZ/2 

V xx 



n 2 v xx dW(x',t) 



(A22) 



T]T]± 



25/4 /\ , {2(x / -r 1 (x)) 2 -r ] 2 V xx ) 
1 dx 



fVx: 



-(i'-,(i)) 2 /, 2 y K 



x e v" '" -"^(aj'.t) (A23) 



Kl + 1 



- 2c 2 (2 (x> (p) + Vi p ) F ra dt 



s/1irV xx rjV± V x i 2 
2V4 ^ t dx' Y^li 2 ^' ~ 71 ^ ~ 



7r 1 / 4 77 3 y / ^X 



x e 



-(x'-r)(x)) 2 /r) 2 V xl 



y3/2 

* XX 

dW{x\t). (A24) 



There is a further simplification that can be made via 
analysis of the diffusion matrix D = BB T ', where B is 
the matrix of innovation terms for the four coupled SDEs 
(A21)-(A24). This matrix has four rows and an infinite 



number of columns (each column represents a different 
noise, and we have a field of noises across x-space). This 
means that D will be a 4 x 4 matrix. Now while the above 
system of SDEs has a unique diffusion matrix, the choice 
of B is not unique. Put another way, this means we are 
free to choose any B that reproduces the correct diffusion 
matrix. In our case, this freedom allows us to reduce the 
number of independent Wiener processes down to two, 
and will also remove the integrals in the above SDEs. 
To begin, let us calculate the D matrix. We have 
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/. . . v x ~ 1/2 (x' - t? ( x ))e- (x '-^) 2 /v 2 v^ 

2 5/4 ^ • • • (V xp /VM 2 ) (x' - 7? (^e-t^W) 1 ^- 

ttVVV^I ••• ri^Vxx 1 ' 2 (2{x' - r) (x)) 2 - ri 2 V xx ) e-( x ' ~^*)) 2 h 2 v x * 

V • • T 1 (K P /KL /2 ) (2(x / - 7? (x)) 2 - 77 2 14.) e-(''^)W v . 



(A25) 



where there are an infinite number of columns, indexed 
by x' . 
matrix: 



D = f3 



where (3 = a/rjrj±. This diffusion matrix can also be 



Prom eq. 


(A25) we can 


calculate 


the diffusion 


/ v 1/2 

j v XX 


V /V 1/2 

v xp 1 V XX 





^ 


V /v 1/2 

v xp/ v XX 


v x %/v x T 














oy-3/2 


OV X p V xx 


V o 





OV X pV XX 










(A26) 



constructed from D 

( 



B'B 

1/4 



IT 



V 3 



Vxx 
/T/3/4 

o 



V o 



where 




VzVxT ^ 

\f~3Vxp I Vxx 



(A27) 



The matrix ( |A27 ) shows that the same system can be 
modelled with considerably simpler innovations terms. 
This simplification leads to the SDEs (31 )-(34) presented 
in Sec. ED 
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